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ABSTRACT 

We  describe  the  development  of  a  two-dimensional  multimaterial  reactive  Eulerian 
hydrocode  to  model  the  detonation  of  condensed  phase  explosives,  and  use  the  code  to 
model  the  projectile  impact  initiation  of  explosives.  The  code  solves  the  Euler 
equations  for  the  conservation  of  mass,  momentum,  and  energy  for  an  inviscid, 
compressible  fluid  using  the  Flux  Corrected  Transport  algorithm.  The  condensed 
phase  materials  are  modelled  using  a  Mie-Gruneisen  equation  of  state,  while  the 
Becker-Kistiakowsky-Wilson  equation  of  state  is  used  to  describe  the  gaseous 
detonation  products.  The  Forest  Fire  reaction  rate  model  is  used  to  describe  the  rate  of 
energy  release  from  the  explosive,  and  a  modified  Young's  algorithm  is  used  to 
maintain  a  sharp  interface  between  different  materials  on  the  computational  mesh.  The 
code  was  used  to  simulate  the  impact  of  a  cylindrical  steel  projectile  against 
Composition  B  explosive,  and  the  threshold  velocity  for  the  onset  of  detonation  was 
found  to  be  dependent  on  the  diameter  of  the  projectile.  Our  computed  results  are  in 
good  agreement  with  experimental  values,  as  well  as  with  results  obtained  from 
Mader's  2DE  reactive  hydrocode.  The  work  described  here  will  provide  the  ADF  with 
an  improved  methodology  to  assess  bullet/fragment  hazards,  and  assist  in  the 
development  of  an  Insensitive  Munitions  capability. 
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Executive  Summary 

The  Australian  Defence  Force  Insensitive  Munitions  (IM)  policy  has  rated 
bullet/ fragment  impact  as  one  of  the  three  major  hazards  to  munitions  in  the 
Australian  defence  context.  A  key  to  reducing  vulnerabihty  is  having  the  capability  to 
predict  the  response  of  candidate  munitions  when  subjected  to  designated  hazards. 
The  work  reported  here  is  concerned  with  understanding  the  processes  and  reactions 
which  control  the  response  of  explosive  fillings  to  various  types  of  projectile  impact. 
This  will  provide  the  ADF  with  an  improved  methodology  to  assess  buUet/ fragment 
hazards  which  can  be  appHed  to  the  assessment  of  existing  munitions,  and  also 
facilitate  the  ADF  "smart  buyer"  capability  for  IM  items. 

The  report  describes  the  development  of  a  two-dimensional  multimaterial  Eulerian 
hydrocode  to  model  the  effects  of  projectile  impact  on  condensed  phase  explosives, 
and  the  effects  of  the  induced  detonation  on  surrounding  materials.  An  important  part 
of  the  study  was  the  use  of  the  code  to  determine  the  threshold  velocity  separating 
detonation  from  non-detonative  events.  The  code  solves  the  Euler  equations  for  the 
conservation  of  mass,  momentum,  and  energy  for  an  inviscid,  compressible  fluid. 
Operator  splitting  is  used  to  reduce  the  complexity  of  the  two-dimensional  calculation, 
and  the  resulting  one-dimensional  equations  are  solved  using  the  Flux  Corrected 
Transport  (FCT)  algorithm  of  Boris  and  Book.  This  has  fourth  order  phase  accuracy, 
and  an  overall  second  order  accuracy  on  uniform  grids.  A  Mie-Gnmeisen  equation  of 
state  is  used  to  describe  the  non-energetic  condensed  phase  materials,  while  the 
Becker-Kistiakowsky-Wilson  equation  of  state  is  used  to  describe  the  reacted 
explosive,  and  a  modified  Mie-Gnmeisen  equation  of  state  is  used  to  describe  the 
unreacted  explosive.  The  rate  of  energy  release  from  the  energetic  material  is  described 
by  the  Forest  Fire  reaction  rate  model. 

A  modified  Yormg's  algorithm  is  used  to  maintain  a  sharp  interface  between  different 
materials  on  the  computational  mesh.  In  this  scheme  an  interface  within  a  cell  is 
approximated  by  a  straight  line.  A  computational  cell  containing  more  than  one  type 
of  material  uses  information  from  each  of  its  eight  neighbouring  cells  to  determine  the 
exact  points  of  intersection  between  the  material  interfaces  and  the  ceU  walls.  This 
removes  the  need  for  the  specialised  comer  treatments  which  are  needed  when  fully 
one-dimensional  interface  tracking  algorithms  are  used. 

We  describe  the  application  of  the  code  to  simulate  the  impact  of  a  cylindrical  steel 
projectile  against  the  miHtary  explosive  Composition  B.  At  relatively  high  projectile 
velocities  the  impact  produces  a  detonation  of  the  explosive,  while  for  lower  velocities 
a  less  violent  release  of  energy  known  as  a  deflagration  occurs.  The  threshold  velocity 
for  the  onset  of  detonation  is  dependent  on  the  diameter  of  the  projectile,  and  we  have 
used  the  multimaterial  code  to  determine  threshold  velocity  as  a  function  of  projectile 
diameter.  Our  computed  results  are  in  good  agreement  with  experimental  values,  as 
well  as  with  results  obtained  from  Mader's  2DE  reactive  hydrocode. 
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1.  Introduction 


Weapons  Systems  Division  is  currently  engaged  in  an  experimental  and  numerical 
modelling  program  to  study  the  response  of  explosive  filled  ordnance  to  projectile 
attack.  The  work  is  being  conducted  in  response  to  the  adoption  of  an  Insensitive 
Munitions  policy  by  the  Australian  Defence  Force,  and  the  identification  of 
bullet/ fragment  impact  as  one  of  the  three  major  hazards  to  munitions  in  the 
Austrahan  defence  context.  The  report  by  Swinton  and  Chick  [1]  describes  the 
development  of  a  projectile  impact  test  designed  to  find  the  critical  impact  velocity  for 
detonation  of  the  target  explosive,  and  its  application  to  bullet  impact  experiments 
conducted  on  Australian  manufactured  Composition  B  explosive.  The  experiments  use 
cylindrical  copper  projectiles,  and  these  are  fired  against  both  bare  and  covered 
Composition  B.  Experiments  have  been  conducted  on  Composition  B  made  from 
boiled  and  milled  RDX  (which  is  the  current  filling  used  in  AushaHan  ordnance),  as 
well  as  on  Composition  B  made  from  recrystalUsed  RDX,  which  is  to  replace  the  boiled 
and  milled  material  when  stocks  run  out.  The  experiments  have  shown  that  the 
Composition  B  made  from  recrystallised  RDX  is  significantly  less  sensitive  to  projectile 
impact  than  that  made  from  boiled  and  milled  RDX,  and  one  of  the  objectives  of  the 
experimental  work  is  to  investigate  the  controlling  parameters  and  mechanisins  of  the 
initiation  process. 

This  report  describes  the  initial  stages  in  the  development  and  appHcation  of  a  two- 
dimensional  (2D)  multimaterial  Eulerian  hydrocode  to  model  the  response  of 
explosive  fillings  to  various  types  of  projectile  impact.  Starkenberg  et  al.  [2]  have 
previously  used  the  Los  Alamos  2DE  code  to  model  projectile  impact  shock  initiation 
of  both  bare  and  covered  Composition  B  charges  using  cylindrical  steel  projectiles. 
They  used  the  2DE  code  to  determine  the  critical  impact  velocity  as  a  function  of 
projectile  diameter  for  bare  explosive,  and  found  that  their  results  were  in  excellent 
agreement  with  the  experimental  results  of  Slade  and  Dewey  [3],  and  the  Jacobs- 
Roslund  empirical  formula  [4]. 

The  2DE  code  is  a  reactive  two-dimensional  Eulerian  hydrocode  which  uses  the  HOM 
equation  of  state  to  model  both  reacted  and  unreacted  explosives,  and  the  Forest  Fire 
reaction  rate  model  to  describe  explosive  decomposition  [5].  The  AMRL  code  under 
development  in  WSD  uses  the  same  equations  as  the  2DE  code  to  describe  the 
explosive  response,  but  the  solution  of  the  transport  equations,  and  the  treatment  of 
material  interfaces,  use  more  advanced  numerical  techniques.  The  AMRL  code  uses 
the  Flux  Corrected  Transport  (FCT)  algorithm  to  solve  the  Euler  equations  [6],  and 
material  interfaces  are  maintained  using  a  Youngs  interface  tracking  algorithm  [7]. 

The  experiments  of  Swinton  and  Chick  have  concentrated  on  the  effect  of  RDX  particle 
size  on  the  critical  velocity  threshold.  The  Forest  Fire  reaction  rate  model  was  not 
designed  to  include  particle  size  effects  however,  and  so  we  have  used  the  AMRL 
code,  in  its  current  state  of  development,  to  model  the  original  experiments  of  Slade 
and  Dewey. 
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In  Section  2  we  outline  the  development  of  the  AMRL  2D  Eulerian  multimaterial 
hydrocode,  and  in  Section  3  we  use  the  code  to  model  the  projectile  impact 
experiments  of  Slade  and  Dewey.  In  Section  4  we  discuss  some  of  the  insights  obtained 
from  the  simulations,  and  describe  the  further  improvements  which  must  be  made  to 
the  code  to  model  the  projectile  impact  experiments  of  Swinton  and  Chick. 


2.  Description  of  Multimaterial  Hydrocode 


2.1  Basic  Equations 

Detonation  phenomena  occur  on  the  microsecond  time  scale,  and  it  is  usual  to  assume 
that  energy  transport  by  heat  conduction,  viscosity,  and  radiation  is  negligibly  small 
compared  with  transport  by  motion.  The  pressures  generated  by  the  detonation  of  a 
condensed  phase  explosive  are  of  the  order  of  10^  atmospheres.  Under  these  conditions 
the  strength  of  any  confining  material  is  completely  negligible,  and  the  material 
responds  hydrodynamically.  In  this  case,  the  appropriate  equations  which  describe  the 
material  response  are  the  Euler  equations,  which  describe  the  conservation  of  mass, 
momentum  and  energy  for  an  inviscid,  compressible  fluid.  These  have  the  following 
form  [8]: 


|^  +  V-(pv)=0 

3t 

(1) 

^^P''^+V-(pvv)=  VP 
dt 

(2) 

dE 

— +  V'(Ev)  =  -V-(Pv) 
dt 

(3) 

In  equations  (1)  -  (3),  p  is  the  density,  v  the  fluid  velocity,  P  the  hydrodynamic 
pressure,  and  E  the  total  energy  per  unit  volume,  which  is  given  by  the  following 
expression 


E  =  pc  +  0.5pv2  (4) 

where  e  is  the  specific  internal  energy.  The  non-reactive  materials  are  described  by  the 
Mie-Gruneisen  equation  of  state,  which  has  the  form  [9]: 

P  =  Ph  +  'y(c-eH)p  (5) 

where  Ph  and  ch  refer  to  the  pressure  and  specific  internal  energy  along  the  shock 
compression  curve,  known  as  the  Hugoniot,  and  which  have  the  form 
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Ph  =  Po  +  poC^oT)/  (1-SoTi)2  (6) 

Ch  =  Co  +  0.5(Ph  +Po)t|/  Po  (7) 

where  r|  is  the  compression,  defined  by  ti  =  1-po/  p,  and  co  is  the  soxmd  speed  in  the 
undisturbed  material.  In  the  above  expressions  the  shock  velocity  Us  and  the  particle 
velocity  Up  are  assumed  to  obey  a  linear  relationship  of  the  form 


tls  —  Co  SoXlp. 


The  energetic  materials  require  different  equations  of  state  to  describe  the  condensed 
phase  material  before  detonation,  and  the  gaseous  detonation  products  after  reaction. 
The  Mie-Gruneisen  EOS  is  used  to  describe  the  unreacted  material,  while  the  BKW 
EOS  is  used  to  describe  the  highly  non-ideal  gaseous  products  after  reaction.  This  has 
the  form  [10]: 


PgVg/(RT)=l  +  x*exp((3x)  (8) 

where  x  =  k/[vg(T+9)“]  and  k  is  the  average  covolume  defined  in  terms  of  the 
individual  covolumes  ki  as  k  =  kSxi  ki,  where  Xi  is  the  mole  fraction  of  species  i.  The  a, 
P,  K,  and  0  are  constants  adjusted  to  reproduce  the  detonation  pressure  and  velocity 
obtained  experimentally. 

The  condensed  explosive  is  converted  into  gaseous  detonation  products  according  to 
the  Forest  Fire  reaction  rate  law  [10], 

^  =  (l-^)exp  EaiP’  (9) 

dt  ^  i=0  2 

where  X  represents  the  fraction  of  reacted  explosive.  The  finite  rate  of  burning 
expressed  by  equation  (9)  yields  a  reaction  zone  of  finite  thickness  in  which  X  varies 
between  0  and  1,  and  which  contains  a  mixture  of  condensed  explosive  and  gaseous 
products.  In  such  regions  we  assume  thermal  and  mechanical  equilibrium  between  the 
two  phases  of  the  explosive  and  use  an  iterative  technique  to  determine  the  specific 
internal  energy  of  each  phase. 


2.2  Numerical  Scheme 

We  have  investigated  two  approaches  to  the  solution  of  the  coupled  equations 
described  above.  In  both  cases  operator  splitting  was  used  to  decouple  the  transport 
stage  from  the  chemical  reaction  stage,  and  then  the  two-dimensional  transport 
equations  were  further  decoupled  into  two  sets  of  one-dimensional  equations.  This 
standard  technique  is  described  in  detail  in  Oran  and  Boris  [6]. 

In  the  first  version  of  the  code  the  density  of  each  material  on  the  grid  was  convected 
independently  by  making  repeated  calls  to  the  LCPFCT  algorithm  [11]  for  each 
material.  The  momentum  flux  and  total  energy  flux  (ie.  specific  internal  energy  plus 
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kinetic  energy)  for  each  material  were  then  added  and  the  combined  fluxes  convected 
using  LCPFCT.  This  approach  is  illustrated  in  the  equations  listed  below,  which  show 
the  sequence  of  equations  to  be  solved  for  the  x-pass  in  a  2D  cartesian  grid  calculation; 
Pi  is  the  density  of  the  i'th  material  on  the  grid,  pu  and  pv  are  the  x  and  y  components 
of  the  TOTAL  momentum  in  each  cell,  and  E  is  the  TOTAL  energy  per  unit  volume  in 
each  cell. 


dPi 

(10) 

dt  dx 

dpu  ^  dpuu  dP 

(11) 

dt  dx  dx 

dpv  1  dpvu  _  Q 

(12) 

dt  dx 

dE  ^  dEu  dPu 
dt  dx  dx 

(13) 

Because  of  the  multimaterial  capability  of  the  code  and  the  diffusive  nature  of  all 
numerical  transport  schemes  it  is  important  to  include  an  interface  tracking  algorithm 
to  maintain  a  sharp  interface  between  different  materials  on  the  grid.  In  the  first 
version  of  the  code  we  used  the  Simple  Line  Interface  Calculation  (SLIC)  scheme  of 
Noh  and  Woodward  [12].  SLIC  is  a  one-dimensional  alternating  direction  method  for 
the  geometric  approximation  of  fluid  interfaces.  The  scheme  constructs  a 
representation  of  the  interface  between  two  materials  from  the  volume  fractions  of  the 
materials  in  the  mixed  cell,  and  by  testing  whether  or  not  the  various  fluids  in  the 
mixed  cell  are  present  or  absent  in  the  zone  just  to  the  left  and  to  the  right  in  the 
coordinate  direction  under  consideration.  The  interface  between  two  different 
materials  is  therefore  represented  by  a  number  of  one-dimensional  components,  each 
of  which  is  composed  entirely  of  straight  lines  either  perpendicular  or  parallel  to  the 
coordinate  direction  under  consideration.  LCPFCT  assumes  that  the  different  materials 
in  a  mixed  cell  are  spread  homogeneously  throughout  the  cell,  and  hence  would 
calculate  an  incorrect  flux  in  a  mixed  cell.  The  SLIC  algorithm  determines  the  correct 
amount  of  material  to  be  advected  into  and  out  of  a  mixed  cell,  and  then  the  LCPFCT 
algorithm  is  modified  by  imposing  a  multiplicative  area  factor  at  the  correct  cell 
boundary  for  each  material.  Further  details  of  this  procedure  can  be  found  in  the 
report  by  Milne  and  Carnegie  [13]. 

If  LCPFCT  and  SLIC  are  applied  to  the  Euler  equations  using  the  above  scheme  then 
problems  are  found  to  occur  in  the  vicinity  of  the  interface  because  the  jump 
conditions  (which  imply  continuity  of  particle  velocity  and  pressure  across  the 
interface)  are  not  enforced.  The  simplest  constraint  which  can  be  applied,  and  is 
consistent  with  the  split  operator  approach  of  FCT,  is  to  ensure  continuity  of  velocity 
across  the  interface.  This  constraint  has  been  added  by  ensuring  that  the  velocities 
either  side  of  a  mixed  cell  are  equal,  with  the  velocity  ahead  of  the  interface  being  set 
to  that  behind  it. 
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A  further  problem  with  the  above  scheme  is  that  the  code  requires  a  unique  pressure 
to  be  deterimned  for  each  mixed  cell,  and  to  do  this  requires  knowing  the  density  and 
specific  internal  energy  of  each  material  in  the  cell.  Whilst  the  above  scheme  convects 
the  density  of  each  material  independently,  only  the  total  energy  in  each  cell  is 
convected,  and  so  a  way  of  dividing  the  total  internal  energy  in  each  cell  between  the 
component  materials  must  be  devised.  To  do  this  we  assumed  pressure  equilibrium 
between  the  materials  in  the  mixed  cell  (this  is  a  common  assumption  in  multimaterial 
Eulerian  codes)  and  used  an  iterative  procedure  to  determine  the  internal  energy  of 
each  component. 

The  above  scheme  was  used  to  simulate  the  one-dimensional  impact  of  a  50  mm  thick 
aluminium  flyer  plate  with  a  velocity  of  1.0  mm/ps  hitting  a  50  mm  thick  slab  of  non¬ 
reactive  PBX-9404  explosive.  The  computational  cell  size  was  1.0  mm,  the  time  step 
was  0.04  ps,  and  the  interface  between  the  aluminium  and  PBX-9404  was  initially 
located  between  cells  50  and  51.  The  pressure  and  particle  velocity  profiles  after  8  ps 
are  shown  in  Figure  1.  These  show  a  shock  moving  forward  into  the  PBX-9404  and 
backwards  into  the  aluminium  flyer.  Simple  shock  impedance  calculations  (Meyers, 
[9])  can  be  used  to  show  that  the  pressure  and  particle  velocity  across  the  interface 
should  have  values  of  0.048  Mbars  and  0.70  mm/ps  respectively,  and  Figure  1  shows 
that  the  simulation  reproduces  these  values  quite  accurately.  Figure  la  however  shows 
that  there  is  a  slight  perturbation  to  the  velocity  profile  at  the  position  of  the  interface 
(which  is  located  partway  through  cell  56),  and  Figure  lb  shows  an  even  greater 
perturbation  to  the  pressure  profile.  While  the  magnitude  of  this  perturbation  is 
relatively  small  (approximately  8%  for  the  pressure),  we  envisage  applications  for  the 
code  where  perturbations  to  the  pressure  amplitude  of  this  size  may  lead  to  larger 
errors,  especially  when  the  code  is  used  to  study  detonation  phenomena. 


Figure  1.  Simulated  (using  SLIC  interface  tracker)  pressure  and  particle  velocity  profiles  after 
8  ps  showing  the  one-dimensional  impact  of  a  50  mm  thick  aluminium  flyer  plate  with  a 
velocity  of  1.0  mmlps  hitting  a  50  mm  thick  slab  of  non-reactive  PBX-9404  explosive. 
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The  source  of  these  disturbances  can  be  traced  to  the  diffuse  treatment  of  the  total 
energy,  and  a  discrepancy  between  the  mass  fluxes  and  interface  velocities  used  in  the 
partial  density  calculations  and  the  momentum  fluxes  used  in  the  principle 
momentum  component  equation.  To  improve  the  performance  of  the  code  it  was 
decided  to  solve  for  the  internal  energy  of  each  of  the  materials,  and  to  replace  the 
SLIC  interface  tracking  algorithm  with  a  scheme  based  on  Youngs  method  [7].  The 
advantage  of  following  the  internal  energy  of  each  material  is  that  an  iterative 
calculation  is  no  longer  needed  in  each  mixed  cell,  and  the  pressure  in  the  cell  is 
determined  from  a  simple  volume  average  of  the  partial  pressures  within  the  cell.  The 
Youngs  interface  tracking  scheme  uses  information  from  all  eight  of  the  mixed  cells 
nearest  neighbours,  and  hence  provides  a  more  accurate  representation  of  the  interface 
than  the  SLIC  scheme.  Additional  changes  to  the  code  included  a  change  in  the  way 
the  principal  component  of  the  momentum  was  calculated.  In  the  new  scheme  the 
momentum  is  calculated  by  using  the  sum  of  the  products  of  the  mass  fluxes 
calculated  during  the  solution  of  the  partial  density  equations,  and  the  cell  interface 
velocities,  also  used  for  the  density  calculations,  thus  ensuring  consistency  in  terms  of 
mass  and  velocity  at  the  cell  boundary.  The  sequence  of  equations  to  be  solved  for  the 
x-pass  in  a  2D  cartesian  grid  calculation  now  has  the  form  shown  below,  where  pi  and 
Ci  are  the  density  and  specific  internal  energy  respectively  of  the  i'th  material  on  the 
grid,  and  u  and  v  are  the  x  and  y  components  of  the  velocity  in  each  cell.  This 
approach  removes  the  need  to  enforce  continuity  of  particle  velocity  across  the 
interface.  A  detailed  description  of  the  new  scheme  can  be  found  in  the  report  by 
Carnegie  [14]. 


^Pi  I  ^  Q 

dt  dx 


(14) 


3u  ^  3uu  _  y 
dt  dx  dx  p 

3v  ^  dvu  _ 
dt  dx  dx 


(15) 

(16) 


de-  ^  de-u  _  p  3u 
dt  dx  dx 


(17) 


The  effect  of  these  changes  can  be  seen  in  Figure  2,  which  shows  pressure  and  particle 
velocity  profiles  for  the  same  ID  simulation  of  a  50  mm  thick  aluminium  flyer  plate 
with  a  velocity  of  1.0  mm/ps  hitting  a  50  mm  thick  slab  of  non-reactive  PBX-9404 
explosive.  The  new  scheme  has  calculated  the  general  features  of  the  flow  quite 
accurately,  although  there  are  differences  between  the  two  schemes.  Figure  2  shows 
that  both  the  particle  velocity  and  pressure  are  continuous  across  the  interface  and  do 
not  display  the  sharp  discontinuities  shown  in  Figure  1.  This  is  quite  a  pleasing  result, 
given  that  the  new  scheme  neither  imposes  a  velocity  constraint  across  the  interface, 
nor  requires  an  iteration  of  the  internal  energy  for  the  materials  in  a  mixed  cell  to 
insure  pressure  equilibrium  between  the  materials.  Figure  2  also  shows  that  the  new 
scheme  provides  a  better  representation  of  the  shock  in  the  aluminium  flyer.  The 
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particle  velocity  plot  in  Figure  2a  shows  a  sharper  shock  profile,  as  does  the  pressure 
profile  in  Figure  2b,  although  there  is  a  slight  perturbation  to  the  pressure  at  the 
location  of  the  front. 


Figure  2.  Simulated  (using  Youngs  interface  tracker)  pressure  and  particle  velocity  profiles 
after  8  ps  showing  the  one-dimensional  impact  of  a  50  mm  thick  aluminium  flyer  plate  with  a 
velocity  of  1.0  mm/ps  hitting  a  50  mm  thick  slab  of  non-reactive  PBX-9404  explosive. 


We  have  also  repeated  this  simulation  using  a  steel  flyer  plate  impacting  non-reactive 
Composition  B.  A  50  mm  thick  steel  flyer  plate  travelling  at  1.0  mm/ps  impacted  a  50 
mm  thick  slab  of  unreactive  Composition  B  explosive.  The  length  of  the  computational 
cell  was  1.0  mm,  the  time  step  was  calculated  from  the  Courant  condition 
(approximately  0.1  ps),  and  the  interface  was  initially  located  between  cells  50  and  51. 
The  results  of  the  calculation  after  8  ps  are  shown  in  Figure  3.  One-dimensional  shock 
impedance  calculations  predict  that  the  pressure  and  particle  velocity  across  the 
interface  should  be  0.057  Mbars  and  0.853  mm/ps  respectively.  Figure  3  shows  that  the 
simulated  particle  velocity  is  approximately  0.84  mm/ps,  while  the  simulated  pressure 
is  approximately  0.057  Mbars.  These  results  are  in  good  agreement  with  the  ID  shock 
impedance  calculations;  the  particle  velocity  is  accurate  to  1%,  while  the  simulated 
pressure  is  no  more  than  5%  too  high.  Figure  3  also  shows  that  both  particle  velocity 
and  pressure  are  continuous  across  the  interface.  Comparison  of  Figure  3  with  Figure  2 
however  shows  a  slightly  greater  perturbation  to  the  shock  profiles  for  the  steel  impact 
problem,  particularly  for  the  pressure  profile.  This  is  caused  by  the  greater  difference 
in  density  between  the  two  impacting  materials  for  the  steel/Composition  B  case  than 
for  the  aluminium/PBX-9404  calculation.  These  perturbations  are  not  fundamental 
however,  and  experimentation  has  shown  that  they  can  be  reduced  by  using  a  finer 
time  step  during  the  first  few  microseconds  of  the  impact. 


7 


DSTO-TR-0325 


Figure  3.  Simulated  (using  Youngs  interface  tracker)  pressure  and  particle  velocity  profiles 
after  8  ps  showing  the  one-dimensional  impact  of  a  50  mm  thick  steel  flyer  plate  with  a  velocity 
of  1.0  mm/ps  hitting  a  50  mm  thick  slab  of  non-reactive  Composition  B  explosive. 


3.  Application  to  Projectile  Impact  Experiments 

Since  the  Forest  Fire  reaction  rate  model  is  unable  to  model  particle  size  effects,  and 
therefore  cannot  be  used  to  simulate  the  Swinton  and  Chick  results,  it  was  decided  to 
test  the  capabilities  of  the  code  in  its  present  stage  of  development  by  modelling  the 
earlier  projectile  impact  experiments  of  Slade  and  Dewey.  These  experiments  used 
steel  cylinders  of  varying  diameters  fired  against  Composition  B  to  determine  the 
critical  impact  velocity  for  initiation  as  a  function  of  the  projectile  diameter.  The 
experiments  were  later  modelled  by  Starkenberg  et  al.  using  the  2DE  code,  and  good 
agreement  was  found  between  the  calculations  and  the  experimental  results. 
Considering  that  the  equations  of  state  and  reaction  rate  law  for  Composition  B  used 
in  our  current  version  of  the  AMRL  code  are  the  same  as  those  used  by  Starkenberg  et 
al.,  we  would  expect  that  simulations  of  projectile  impact  using  our  code  should  be  in 
good  agreement  with  the  results  found  by  the  above  authors. 

3.1  Geometry  and  Computational  Considerations 

The  projectile  impact  calculations  were  performed  using  2D  cylindrical  geometry.  The 
basic  grid  is  illustrated  in  Figure  4.  This  shows  the  cylindrical  steel  projectile,  the 
cylindrical  Composition  B  explosive  charge,  and  the  surrounding  air.  In  all  the 
simulations  the  projectile  was  initially  just  in  contact  with  the  explosive  and  had  a 
positive  initial  velocity,  indicating  that  the  direction  of  motion  was  towards  the 
explosive.  The  initial  pressure  of  the  air,  projectile  and  high  explosive  was  one 
atmosphere  (ie.  101325Pa),  and  this  determines  the  state  of  the  respective  materials  via 
the  appropriate  equation  of  state. 
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Figure  4.  Schematic  illustration  of  the  2D  cylindrical  geometry  used  in  the  calculations 
showing  the  cylindrical  steel  projectile,  the  cylindrical  Composition  B  explosive  charge,  and  the 
surrounding  air. 


3.2  Steel  projectile  /  Bare  Composition  B  impact  results 

The  critical  velocity  (Vc)  is  defined  as  the  velocity  of  the  projectile  above  which 
projectile  impact  results  in  shock  initiation,  and  below  which  the  impact  results  in  a 
lower  order  event  such  as  a  burn  or  deflagration.  Experimentally  it  has  been  found 
that  the  critical  velocity  increases  as  the  diameter  of  the  projectile  is  reduced.  In  the 
calculations  reported  here  the  speed  of  the  projectile  has  been  varied  to  determine  the 
critical  velocity  for  each  projectile  diameter.  In  this  process,  it  is  important  to  carefully 
determine  the  type  of  event  caused  by  the  projectile  impact.  Starkenberg  et  al.  [2]  used 
contours  of  the  explosive  mass  fraction  to  determine  the  outcome  of  projectile  impact 
simulations,  with  detonation  being  recognised  by  the  close  spacing  of  the  contour 
lines.  In  this  context,  mass  fraction  refers  to  the  amount  of  unreacted  explosive,  thus  a 
mass  fraction  of  1  corresponds  to  wholly  xmreacted  explosive,  and  a  mass  fraction  of 
zero  represents  complete  reaction  of  the  explosive.  In  this  work,  colour  contour  plots  of 
the  pressure  and  the  speed  of  the  reaction  front  are  also  used  to  determine  the  outcome 
of  the  projectile  impact  calculation. 

Figure  5  shows  the  results  from  the  impact  of  a  6.75  mm  diameter  steel  projectile  on 
bare  Composition  B.  In  this  example  the  projectile  velocity  is  1.24  mm/|is  and  the 
impact  results  in  a  detonation.  This  is  evident  from  the  close  spacing  of  the  black 
contour  lines  of  the  explosive  mass  fraction,  and  by  the  high  pressure  at  the  shock 
front  shown  by  the  colour  pressure  contours  (red  represents  the  highest  pressure). 
Figure  6  shows  contour  plots  from  an  identical  run,  except  with  a  projectile  velocity  of 
1.12  mm/ps,  and  this  clearly  shows  a  failure.  In  this  case  the  black  contour  lines  of  the 
explosive  mass  fraction  are  widely  spaced  and  the  pressure  at  the  shock  front  is 
relatively  low. 


9 


DSTO-TR-0325 


Time=  1.0938/2,3  Cycle—  451 


Time=  1.3846/2,3  Cycle=  601 


.M-04  .l«-«  jm-m  .M»-«  .»-« 


tm-m  .m-*i  M»-m  ja«-«  .i»^  -m-tt  .m-n 


JOkA  JB-41  .IMMB  .U>t40  i«t<W 


Time=  2.1938/2,3  Cycle—  1051 


Time=  2.6515/2,3  Cycle—  1501 


Figure  5.  Simulated  results  from  the  impact  of  a  6.75  mm  diameter  steel  projectile  travelling  at 
1.24  mm/ms  on  bare  Composition  B. 
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Time=  1.1886yLts  Cycle=  451  Time=  1.5373yU.s  Cycles  601 


Time=  2.5071/2S  Cycle=  1051  Time=  2.7575/2S  Cycle=  1501 


Figure  6:  Simulated  results  from  the  impact  of  a  6.75  mm  diameter  steel  projectile  travelling  at 
1.12  mm/jus  on  bare  Composition  B. 
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Figure  7.  Plot  of  shock  front  position  as  a  function  of  time  for  the  6.75  mm  projectile  impacting 
the  bare  explosive  at  2.24  mm/ps.  The  dots  denote  the  computed  positioUj  while  the  solid  line  is 
a  line  of  best  fit  calculated  from  the  last  third  of  the  points,  indicating  that  the  speed  of  the 
shock  front  increases  towards  the  end  of  the  run. 


The  velocity  of  the  impact  shock  front  also  gives  an  indication  of  the  type  of  outcome 
produced  by  the  projectile  impact.  The  plot  of  shock  front  position  as  a  function  of 
time  for  the  6.75  irun  projectile  impacting  the  bare  explosive  at  1.24  mm/  ps  is  shown  in 
figure  7.  The  slope  of  this  plot  represents  the  speed  of  the  shock  front  and  shows  that 
as  reaction  proceeds  the  speed  of  the  shock  front  increases  and  reaches  a  value  of 
about  7.1  mm/ps.  On  the  other  hand,  the  corresponding  plot  for  the  6.75  mm  projectile 
travelling  at  1.12  mm/ps,  shown  in  Figure  8,  exhibits  the  reverse  behaviour.  In  this 
case  the  velocity  of  the  shock  front  decreases  and  the  final  velocity  is  approximately  4.2 
mm/ps.  Although  the  detonation  velocity  of  Composition  B  is  8.0  mm/ ps,  the  value  of 
7.1  mm/ps  is  not  the  final,  steady-state  velocity.  If  the  calculation  was  allowed  to 
proceed  further,  the  reaction  front  velocity  would  increase  towards  the  expected  value. 
Nonetheless,  plots  of  shock  front  position  as  a  function  of  time  help  to  determine  the 
outcome  of  the  projectile  impact.  If  the  shock  front  velocity  increases  with  time,  the 
reaction  is  accelerating  and  this  indicates  a  detonation,  whereas  a  decelerating  shock 
front  is  indicative  of  a  failure. 

Table  1  summarises  the  important  parameters  of  the  computations  performed  in  our 
study  of  steel  projectile  impact  on  bare  Composition  B.  In  this  table,  d  is  the  projectile 
diameter,  dx  and  dy  are  the  grid  sizes  in  the  x  and  y  directions,  ncx  and  ncy  are  the 
number  of  cells  in  the  x  and  y  directions.  Cycles  is  the  number  of  time  steps  performed 
in  the  simulation,  and  Elapsed  Time  is  the  cumulative  total  of  dt,  the  time  step,  at  the 
end  of  the  nm.  The  outcome  of  the  particular  projectile  impact  is  deduced  by 
examining  the  mass  fraction  contours,  colour  contour  plots  of  the  pressure  profiles, 
and  plots  of  the  shock  front  position  as  a  hmction  of  time. 
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Figure  8.  Plot  of  shock  front  position  as  a  function  of  time  for  the  6.75  mm  projectile  impacting 
the  bare  explosive  at  1.12  mm/fis.  The  dots  denote  the  computed  position,  while  the  solid  line  is 
a  line  of  best  fit  calculated  from  the  last  third  of  the  points,  indicating  that  the  speed  of  the 
shock  front  decreases  towards  the  end  of  the  run. 


The  critical  velocity  is  known  to  be  an  inverse  function  of  the  projectile  diameter,  ie., 
Vc=  kd'^' .  Starkenberg  et  al.  have  used  n  =  0.5  to  plot  their  results  on  the  projectile 
impact  of  Composition  B,  but  n  =  1  appears  to  give  a  better  fit  to  the  experimental  data, 
as  demonstrated  by  Liddiard  and  Rosland  [15].  Hence  in  this  report  we  use  n  =  1  when 
comparing  our  results  with  previously  published  work. 

Figure  9  shows  a  plot  of  Vc  versus  1/d  for  our  results,  the  experimental  results  of  Slade 
and  Dewey  [3],  and  the  2DE  computations  of  Starkenberg  et  al.  [2].  The  2DE  code  uses 
the  HOM  EOS  and  the  Forest  Fire  reaction  rate  model  to  describe  the  explosive 
response,  but  the  treatment  of  both  the  material  transport  and  material  interfaces  is 
entirely  different  from  the  methods  described  here.  2DE  uses  a  finite  difference  scheme 
based  on  the  work  of  Gentry  et  al.  [19]  to  solve  the  transport  equations.  This  is  a  low 
order  scheme  which  relies  on  the  use  of  an  artificial  viscosity  to  stabilize  the  shock,  and 
hence  is  less  accurate  than  the  FCT  scheme  employed  here.  The  treatment  of  mixed 
cells  in  2DE  is  based  on  the  donor-acceptor  method  developed  by  Johnson  [20].  This 
has  limited  resolution,  and  is  tmable  to  locate  the  position  of  the  interface  within  the 
mixed  ceU.  The  Yoimgs  interface  tracking  algorithm  used  in  the  present  code  has 
greater  resolution,  and  is  able  to  draw  an  accurate  representation  of  the  interface 
within  each  cell. 
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Table  1:  Steel  projectile  impact  on  bare  Composition  B  -  summary  of  runs. 


„d  . 

(mm) 

dx/dy 

(mm) 

ncx 

ncy 

Cycles 

Elapsed  Time 
(ps) 

Projectile  velocity 
(mm/  ps) 

Outcome 

4.0 

0.133 

200 

60 

4500 

2.481 

1.66 

Fail 

4.0 

0.133 

300 

60 

6000 

3.120 

1.68 

Det. 

5.0 

0.167 

200 

60 

4500 

2.929 

1.38 

Fail 

5.0 

0.167 

200 

60 

6000 

3.037 

1.40 

Det. 

6.75 

0.225 

160 

60 

3000 

4.450 

1.14 

Fail 

6.75 

0.225 

160 

60 

3000 

3.981 

1.16 

Det. 

8.0 

0.267 

120 

60 

1500 

3.299 

1.04 

Fail 

8.0 

0.267 

120 

60 

1500 

3.458 

1.06 

Det. 

10.0 

0.333 

200 

60 

3000 

4.575 

0.92 

Fail 

10.0 

0.333 

120 

60 

1500 

4.685 

0.94 

Det. 

12.0 

0.4 

200 

60 

3000 

6.948 

0.82 

Fail 

12.0 

0.4 

200 

60 

3000 

6.925 

0.84 

Det 

15.0 

0.5 

200 

60 

3000 

7.007 

0.72 

Fail 

15.0 

0.5 

120 

60 

1500 

6.939 

0.76 

Det. 

18.0 

0.6 

300 

60 

1500 

7.955 

0.67 

Fail 

18.0 

0.6 

300 

60 

1500 

8.124 

0.68 

Det. 

Our  calculated  results  are  in  close  agreement  with  the  calculations  of  Starkenberg  et 
al.,  although  there  appears  to  be  a  consistent  discrepancy  at  larger  projectile  diameters, 
where  the  Starkenberg  et  al.  results  indicate  a  slightly  lower  critical  velocity  in 
comparison  to  our  results.  The  magnitude  of  the  difference  is  not  great  however,  and 
we  believe  this  may  be  attributed  to  the  different  methods  used  to  distinguish  between 
detonation  and  failure  by  the  various  authors,  rather  than  to  the  differences  between 
our  code  and  the  2DE  code.  It  should  be  noted  that  in  our  work  the  difference  in  the 
velocity  of  the  projectile  that  produces  a  fail  and  a  detonation  is  typically  0.02  mm/ ps. 
Starkenberg  et  al.  do  not  specify  the  change  in  projectile  velocity  between  a  fail  and  a 
detonation.  The  discrepancy  between  either  of  the  computed  critical  velocities  and  the 
experimental  results  is  larger  than  the  difference  between  the  two  sets  of  computed 
values.  The  calculated  results  consistently  predict  a  lower  critical  velocity  than  the 
experimental  results.  This  can  be  attributed  to  deficiencies  in  the  reaction  rate  model. 
Forest  Fire,  used  for  the  explosive. 
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Figure  9.  Plot  ofVc  versus  1/d  for  our  results,  the  experimental  results  of  Slade  and  Dewey, 
and  the  2DE  computations  of  Starkenberg  et  al. 


4.  Discussion  and  Conclusion 


The  simulation  results  shown  in  Figures  5  and  6  are  examples  of  either  a  clear  failure, 
or  a  clear  detonation.  However,  during  the  course  of  the  calculations  it  became 
apparent  that  the  result  was  not  always  immediately  clear,  and  it  became  necessary  to 
perform  longer  runs  to  clarify  the  outcome.  Figure  10  shows  another  example  where  a 
clear  failure  is  immediately  apparent.  In  this  case  the  projectile  had  a  diameter  of  12.0 
mm  and  an  impact  velocity  of  0.76  mm/ ps.  In  each  of  the  four  plots  the  pressure  levels 
remain  low  and  the  mass  fraction  contour  lines  always  remain  well  spaced,  and  none 
of  the  plots  show  any  indication  that  a  detonation  will  develop.  The  sequence  shown 
in  Figure  11  also  represents  a  failure,  however  in  this  case  it  initially  looks  like  a 
detonation.  The  projectile  had  a  diameter  of  4.0  mm,  and  an  impact  velocity  of  1.60 
mm/ ps.  At  0.3203  ps  there  is  a  region  of  high  pressure  just  behind  the  reaction  front 
and  the  mass  fraction  contour  lines  are  closely  spaced,  but  at  0.9295  ps  the  mass 
fraction  contour  lines  have  spread  out  indicating  that  the  reaction  is  slowing  down.  By 
1.8336  ps  the  contour  lines  become  even  further  spaced  and  the  pressure  has  reduced 
significantly,  as  indicated  by  the  colour  contours.  This  type  of  behaviour,  where  the 
detonation  dies  away,  is  observed  for  the  smaller  diameter  projectiles. 

In  other  cases  the  initial  stages  may  look  like  a  failure,  but  ultimately  a  detonation 
results.  Such  a  sequence  is  shown  in  Figure  12,  where  the  diameter  of  the  projectile  is 
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12.00  mm  and  the  impact  velocity  is  0.84  mm/ps.  At  4.4090  ps  the  mass  fraction 
contour  lines  are  widely  spaced  and  the  colour  contours  indicate  relatively  low 
pressures.  However/by  6.0114  ps  a  strong  reaction  front  has  developed  with  closely 
spaced  mass  fraction  contour  lines  and  high  pressure  behind  the  front.  This  buildup 
continues  and  by  6.5665  ps  a  detonation  has  clearly  formed. 

An  example  of  a  prompt  detonation  sequence  is  shown  in  Figure  13.  In  this  case  the 
projectile  diameter  was  18.0  mm  and  the  impact  velocity  was  0.82  mm/ps.  High 
pressures  and  closely  spaced  mass  fraction  contour  lines  are  apparent  soon  after  the 
projectile  impact.  After  this  immediate  initiation,  the  detonation  continues  to 
propagate.  This  type  of  prompt  initiation  is  typical  for  the  larger  diameter  projectiles. 

The  results  presented  above  have  shown  the  capabilities  of  the  code  in  its  present  state 
of  development,  and  our  ability  to  accurately  predict  the  critical  velocity  for  projectile 
impact  on  Composition  B  when  the  Forest  Fire  coefficients  appropriate  to  the 
particular  composition  under  study  are  used.  It  is  impractical  however  to  develop  the 
correct  Forest  Fire  coefficients  for  each  of  the  compositions  studied  by  Swinton  and 
Chick,  and  so  an  alternative  procedure  will  have  to  be  developed  to  model  these 
results.  Because  the  experiments  show  such  a  strong  dependence  on  RDX  particle  size 
and/or  morphology,  any  reaction  rate  scheme  used  to  model  the  results  will  have  to 
include  an  explicit  dependence  on  particle  size.  Appropriate  reaction  rate  schemes 
currently  do  not  exist;  however,  by  combining  the  Lee  and  Tarver  Ignition  and  Growth 
reaction  rate  model  [16]  with  the  work  of  Kim  [17]  on  particle  size  dependent  reaction 
rates  in  shocked  explosives,  and  the  work  of  Murphy  et  al.  on  Composition  B  shock 
initiation  [18],  we  anticipate  that  some  insight  into  the  experimental  results  might  be 
obtained.  Such  work  is  proceeding. 
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Figure  10:  Simulated  results  from  the  impact  of  a  12.0  mm  diameter  steel  projectile  travelling 
at  0.76  mm/ps  on  hare  Composition  B.  (clear  failure  is  immediately  apparent) 
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Figure  11:  Simulated  results  from  the  impact  of  a  4.0  mm  diameter  steel  projectile  travelling  at 
1.60  mm/jjs  on  hare  Composition  B. 
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Figure  12:  Simulated  results  from  the  impact  of  a  12.0  mm  diameter  steel  projectile  travelling 
at  0.84  mm/ps  on  bare  Composition  B.  (initially  looks  like  a  fail,  but  ultimately  a  detonation 
results). 
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Figure  13:  Simulated  results  from  the  impact  of  a  18.0  mm  diameter  steel  projectile  travelling 
at  0.82  mm/ps  on  hare  Composition  B.  (example  of  a  prompt  detonation  sequence). 
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